NOTICE 


THIS DOCUMENT HAS BEEN REPRODUCED FROM 
MICROFICHE. ALTHOUGH IT IS RECOGNIZED THAT 
CERTAIN PORTIONS ARE ILLEGIBLE, IT IS BEING RELEASED 
IN THE INTEREST OF MAKING AVAILABLE AS MUCH 
INFORMATION AS POSSIBLE 



STRUCTURE AND SCREENING IN MOLECULAR AND METALLIC 
HYDROGEN AT HIGH PRESSURE* 


(NASA-CB-164467) STBOCTOBE AMD SCBEENING IN 
MOLECULAB AND ilEIALLIC HYDBOGl'N AT HIGH 
PBESSUBE (Cornell Oniv., Ithaco;, N. Y.) 

32 p HC A03/MF A0 1 CSCL 20H 

G3/72 


D. M. Wood’^ and N. W. Ashcroft 
Laboratory of Atomic and Solid State Physics 
Cornell University, Ithaca, New York 14853 



4 ' 


N81-25780 

(Jnclas 

27861 





April 1981 
MSC Report #4460 


2 


Abstract 


A variational wavef unction proposed by Abrikosov is used to express 
the (spin-restricted ) Hartree-Fock energy as reciprocal lattice suns for 
static lattice FCC monatomic hydrogen and diatomic Pa3 molecular hydrogen. 

In the monato;:t,ic phase the hydrogenic orbital range is found to closely 
parallel the inverse Thomas- Fermi wavevector; the corresponding energy E 
has a minimLn of -0.929 Ryd /electron at r^ = 1.67. For the diatomic phase 
E(r^) is similar, but the constituent energies, screening, and bond length 
all reflect a qualitative change in the nature of the solid at r^ = 2.8. 

This change is interpreted in terms of a transition from protons as struc- 
tural units (at high density) to weakly interacting molecules (at low den- 
sity). Insensitivity of the total energy to a rapid fall in the bond length 
suggests association with the rotational transition, where the rigid mole- 
cular orientations characteristic of high pressures disappear and the mole- 
cules rotate freely at low pressure. The importance of phonons (neglected 
here) in a correct treatment of the total energy is emphasized, and the 
possible connection between the rotational transition and metallization 
of the diatomic phase is discussed. It is concluded that methods which 
sphericalize the Wigner-Seitz cell may overlook important structural proper- 
ties (to which che total energy is relatively insensitive) for the diatomic 
phase. 
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I : Introduction 

Hydrogen continues to attract attention among condensed matter physicists 
because it has the potential to be the simplest of all possible physically 
realizable metals. Because of its low nuclear mass and the metallic nature 
of the high density condensed phase, metallic hydrogen is expected to exhibit 
interesting electronic ordered states, for example, high temperature super- 
conductivity. ’ Indeed, precisely because dense hydroger' is a quantum 
system--one whose ze: o point energy is significant in comparison to the 
binding energy of the condensed phase itself--it may not be a solid at all, 
but rather a metallic liquid ’ (and possibly the first known liquid super- 
conductor ), with associated interesting two-component Fermi liquid effects 
at low temperatures . For our purposes we will assume the monatomic phase 
to be crystalline, in fact, face centered cubic (in accordance with the 
work of Straus and Ashcroft^). This assumption is reasonable since the 
energies associated with ordering are exceedingly small compared with the 
structural energies that we shall encounter. 

The molecular solid observed for low temperatures and pressures exhibits 
interesting properties of its own, again associated with the low mass of 

O 

the hydrogen atom. At low pressures the remarkable sphericity of para- 
hydrogen molecules and their small moment of inertia (a consequence of the 
compactness of a molecule) permits them to rotate freely. For p-ca-hydrogen 
at higher pressures, where the enhanced electric quadrupole-quadrupole inter- 
action between molecules locks them into orientational order, there are 
elementary excitations in the form of librons, i.ei, quantum mechanical 
"zero point motion" associated with the (not always) small amplitude orienta- 
tional oscillations of the hydrogen dunbbells about their crystalline 
orientations. (The characteristic temperature scale for these excitations 
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is somewhat lower for ortho-hydrogen.) In what follows we shall assune that 
the structure of ortho- hydrogen at low temperatures and pressures (the 
space group Pa3) is appropriate for both ortho- and para-hydrogen at high 
pressures . 

The process of metallization of hydrogen may occur in two distinct 
ways: (i) As in the case of iodine, another covalently bonded diatonHc 

solid, the application of pressure may continuously reduce the overall gap 
between valence and conduction bands to the extent that, without change in 
structure band overlap occurs, thereby giving rise to conduction. The basic 
diatomic order remains, (ii) The hydrogen molecules may be directly 
pressure-ionized and dissociated, the associated change in structure giving 
rise to a monatomic metallic phase. Part of the work to be described here 
is directed toward an understanding of the manner in which this dissociation 
occurs . 

Many past estimates of the equations of state and the transition pres- 
sure from diatomic to monatomic hydrogen have relied on different calculational 
techniques for the two phases. Often, for example, the energy of the molecu- 
lar phase has been taken from a superposition of pair potentials, while 
most recent treatments of metallic hydrogen have used perturbation theory 
about the uniform interacting electron gas.^^ There are, however, funda- 
mental reasons for believing that a pair potential approach will lose validity 
at high pressures (when the energy of molecular interaction becomes compare- 

in ^ n 

ble to the lattice binding energy), requiring a scheme better suited 
to the delocalized nature of the resulting electron wavefunctions . 

We shall adopt below an "exact exchange" Hartree-Fock description of 
the condensed phases of hydrogen. There are, of course, more realistic 
formalisms which are on a firmer footing as regards the i ncorporation of 
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many-body (i,e., correlation) effects, chief among these the density-functional 
schemeJ^ As usually implemented, these procedures include in an approximate 
local exchange-correlation potential the physical effects of electron indis- 
tingui shabi 1 i ty. All such local approximations are enormously simpler compu- 
tationally than the corresponding Hartree-Fock calculations, with their 
characteristic non-local exchange terms. 

However, as yet no density-functional ''alculation has been applied to a 
system of full thcee-dimensional crystalline symmetry, i.e., it has been cus- 
tomary to replace the Wigner-Seitz cell of the crystal with a sphere of equal 

volume, and the periodic boundary conditions of the crystal with the Wigner- 
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Seitz boundary condition. Although in terms of total energy this procedure 
typically gives rise to proportional shifts expected to be of the order of the 
difference between Madelung energies of different cubic structures (0(10"^)), a 
great deal of structural information is relinquished in such a description. This 
is of relatively little consequence for monatomic phases, but for molecular 
hydrogen at high pressures, as has been discussed by Chakravarty et al.,^^ 
phonon energies are quite comparable both to the static energy di fferences 
between different crystal structures and to the solid binding energy per elec- 
tron itself. It is therefore of interest to examine calculational schemes in 
which such small structural differences are retained. While we shall com- 
pletely neglect the effect of lattice vibrations below, we will retain a 
complete description of the diatomic crystalline structure, with the hope of 
understanding the behavior of the structural and screening properties of the 
diatomic phase at high pressures. 

Hartree-Fock (HF) theory for solids, as mentioned above, has been com- 
paratively rarely implemented. ^ It is useful to keep in mind a number 
of limitations of the Hartree-Fock approach: (i) The correlation energy 

(by definition) is omitted. Insofar as we shall be concerned with structural 
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details, and the density dependence of the correlation energy is believed 
i:o be weak, we shall not be concerned about this problem^ (ii) Hartree- 
Fock calculations for the band structures of metals typically give rise 
to bandwidths (i.e., the energy difference between the bottom of the valence 
band and the Fermi level) excessively large in comparison to either competing 
calculational methods or to experiment. For a total structural energy cal- 
culation, such as we perform below, however, the energy per electron generally 
is not unsatisfactory, (iii) Spin-restricted Hartree-Fock, in which each 
electron level has twofold spin-degeneracy, leads to physical curiosities 
in the extreme low-density limit cf a crystal. Instead of reproducing pro- 
perties of an isolated atom, spin-restricted H-F typically gives rise to 
neutral or even partially charged pseudo-atoms which, even if neutral, attract 
each other at long distances.' As noted by Stanton, this is essentially 
because the entity associated with a pseudo-atom is a fraction of an electron 
pair, not a single electron: it is not surprising that the high on-site 

correlation energy (the repulsive energy associated with an electron pair 
on a pseudo-atom) makes the binding energy of a pseudo-atom anomalously 
small. We shall return to this point below in a discussion of the results 
for metallic (monatomic) hydrogen. 


1 1 : Abrikosov Variational Wavef unction 

1 9 

We will adopt a wavefunction proposed by Abrikosov as a reasonable 
choice for a unified description of both molecular and metallic hydrogen 
at high pressure. It is 



it-r 


I cf)(r-’^-'S) 


( 1 ) 
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wi th 


<j)(r) = exp(-Ar) 


( 2 ) 


where we have chosen a (lattice + basis) description of the crystal structure 
of interest, with the underlying Bravais lattice and {t} the set of basis 
vectors within a conventional unit cell. The subscript 1 < is intended to indi- 
cate the delocalized nature of this wavefunction: it is explicitly of Bloch 

form. By the use of a variational parameter A, we have allowed for the possi- 
bility that in a condensed environment the shape of a hydrogen wavefunction may 
differ from that of an isolated atom. In the molecular phase the protons are 
grouped in pairs centered at the sites of the Bravais lattice and both A and 
the molecular inter-proton spacing 2D (which enters straightforwardly in the 
specification of the Pa3 lattice) will be such variational parameters. 

A few remarks about the form of (1) are appropriate. In the standard 
tight binding method one constructs tne wavefunction by superposing linea*,* 
combinations of localized atomic wavefunctions , i.e., 

(3) 

where 

" I (4) 

n 

and the 'j^p(r) satisfy the Schrbdinger equation for an isolated atom; one 
hopes that relatively few n are needed in (4). Our choice (1) differs from 
(3) in that (i ) we are using only one atomic orbital, (ii) the choice of 
phase is si te-i ndepende.nt; this will mean important simplifications later 
in the evaluation of, for example, the kinetic energy. 
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Point (i) is not a serious approximation if one is specifically inter- 
ested in the electronic ground state. The choice of phase (point (ii)) 
allows for delocalization of conduction electrons while, preserving the essen- 
tial atomic-like nature of the wavefunctions near the proton sites. We 
will see that this choice will mean that free electron behavior will emerge 
naturally in the limit of large densities. In connection with the choice 
of phase it is also useful to observe that for the molecular solid the scale 
of significant variation of the wavefunction (1) is set by the molecular 
inter-proton spacing 2D, while the scale of k is such that k kp where 
kp is the Fermi wavevector. For low pressures 2kpD n. 0.15 « 1 so that 
the chosen variational function will tend to preferenti al ly deposit charge 
between the protons in the same molecule, i.e., there is relatively little 
of the "anti bonding" orbital mixed in, and (1) should be satisfactory for 
a variational calculation. 

p 

In particular, it should be noted that since li|^|^(r,)l is independent 
of the Bloch wavevector we are restricted to the description of a system 
with a spherical Fermi surface. It is known from the structural expansion 
method, for example, that effects originating from departures from sphericity 
of the Fermi surface first occur in fourth order in perturbation theory 
about the homogeneous interacting electron gas, the perturbation being 
the electron-ion interaction.^^ This therefore suggests that, at least 
at high pressures, the assumed absence of a dependence of the energy on 
the direction of twill also be a satisfactory approximation. For possible 
cubic phases of (monatomic) metallic hydrogen this assumption is quite well 
justified; for the insulating molecular phase, where one might imagine 

proceeding via a tight-binding calculation, this approximation is probably 
poor at quite low density. 
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III: Hartree-Fock Constituent Energies 

As noted, in what follows we shall neglect the kinetic energy of the 

ions (i.e., phonons) within the fraiticwork of the Born-Oppenhoimer (adiabatic) 

. . 24 

approximation. 

For our purposes a description in momentum space is convenient. With 
Z = 1 and wavefunctions (1) our Hartree-Fock states are specified by the 
Block wavevector and the spin. For a crystal of N cells and volune Q the 
wavef unction becomes 

<5) 

G 


and the crystal potential (pure Coulomb for hydrogen) is 


V(r) = 






( 6 ) 


where nj^ is the number of sites in the basis, and for reciprocal lattice 
vector s is the normalized basis structure factor 


(7) 

The constituent Hartree-Fock energies are kinetic (<K.E.>), electron- 
proton (<V^ >), electron-electron (<E^ ^>) and exchange (<E„„>). Per 

c“P c“C cX 

electron, these are, respectively, 
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is the noma! i zed Gth Fourier component of the electron density, Cp is the 
Fermi energy and kp the Fermi wavevector, and - Nnj^. Here 


sL " /n\s 1 j dq 


(Err) 


4'n’ 
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(13) 


with both integrations extending over a spherical Fermi surface. Explicitly 


f^(x) = - Ex) £n 


2+X 


2-X 




} . 


(14) 
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This function was first evaluated by Abrikosov. 

The total energy also includes the Madelung energy i.e., the energy 

of a lattice of point positive ions (embedded in a uniform background qf 
compensating negative charge) interacting via the^ Coulomb potential. The 
total energy per electron (neglecting phonons) may thus be written 


n 



1_ 3 





(15) 


where E ,(J/0), the electronic energy, 1s the sun of terms (7)-(11) except 
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that m have rewitten the G » 0 contributions to the kinetic and exchange 
energies in terms of the familiar dimensionless density parameter r^ defined by 


4ir 

3 



(16) 


where n is the electron nunber density and ag is the Bohr radius. We then 
recognize the terms in brackets above as the kinetic and exchange energies 
of a uniform interacting electron gas; the ^ = 0 Hartree-Fock energy could 
be supplemented by an estimate of the correlation energy, e.g. , the 
Nozi$res-Pines interpolation formula.^® No =5 0 terms appear in the electro- 
static energies because the system as a whole is neutral. 

For purposes of computation it is convenient to back-transform into 
real space the simple k-space sums appearing in the kinetic energy. This 
is expected to be beneficial since our orbital decays exponentially 
for large distance, so the real space sums we generate should converge well, 
at least at low densities where the system is most non-uniform* We find 


f|2 2 


1 +1- I I + X - X^/3) 
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( 17 ) 
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where X - \|R + b - b‘|. We find, furthennore, 


^ ^nor |i [z^ + (g+g')^3^ Cz^ + g'^]^ 


s*(l‘)s(G+Q') 
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where the nontalizing constant is 
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Note that for 6 = ~ (i,j,k), g (i,j,k). Taken togecher the total energy 
for hydrogen is 


c =, c: + rp + nc 
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( 21 ) 


We have used the general form of the Madelung energy 


^Mad = 


( 22 ) 


In the spirit of the structural expansion mentioned above, we may now 
identify the "band structure energy." It is that portion of the total energy 
due to nufi- uniformity of the electronic system, To a11 orders in the 
electron^ion interaction in this model, it is the term in ( ) above. 

IV: Computational Details 

Equation (20) is the basis for a Hartree-Fock calculation of the zero 
temperature equation of state of both monatomic and molecular (diatomic) 
hydrogen, results of which are presented below. Inspection of Eqs . (20) and 
(18) indicate that the inhomogeneous (6 / 0) contributions to the energy 
involve a sixfold sunmation over the reciprocal lattice indices (i,j,k) and 
(i'J'sk'). It is thus imperative to exploit such symmetries as may be 
present to facilitate evaluation of the sums. Further details are given 
in Appendix A. 

A unit cell of the Pn3 structure, believed to be appropriate to solid 
diatotfiic hydrogen at low temperatures and high pressures for both the ortho 
and para phases is shown in Figure 1. It is most simply described as 
a simple cubic conventional cell (of side a) with an eight-point basis, 
with ions at the points 

T.z 

■S, 4 = a[(0,4,4) ± 

(23) 

bg^g = a[(i,0,i) ± r(l,1,-l)3 

b 7^8 = a[(i,i,0) ± r(-l.l,l)] 

where r = D/(/3a) and 2D is the "bond length" of a given molecule. (The 

structure is FCC with molecules at each site, each oriented along a 

different body diagonal.) This is a highly uns.ymmetrical structure-, its 
normalized basis structure factor for a simple cubic reciprocal lattice 
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vector G = 1^ (i ,j,k) is 


S(i ,j ,k) = 2[cosa(i + j + k) + (-1 )^'^'^cosa(-i + j + k) 


+ (-1 )'^'^^cosa(i - j + k) + (-1 )''^'^cosa(i + j - k)] 


i+k. 


( 24 ) 


where a = . For D/a = /5/4 the ions are equally spaced along the body 

diagonals, and the structure becomes simple cubic with a lattice constant 
a/2. The Madelung energy per electron for this structure may be computed 
using the usual Ewald expression 


F sLUl 
M ‘ 2=0 >-s 


3 I l^|s(^)|2e-sV^l 
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ill 




(25) 


where P is the free parameter resulting from the Ewald process (i.e., Ej^ 
should be numerically independent of P for P of order 1), erfc is the com- 
plementary error function, g = Ga/(27r), and the prime indicates omission 
of the ^ = 0, S = S' term, i.e.. 


I 



Z I + I I 


R/0 b/b' R=0 b/b' 



The Madelung constant for the Pa3 structure, as a function of the 
parameter D/a, is displayed in Figure 2; also shown for comparison is the 
Madelung constant for a "rotationally averaged Pa3" structure, discussed 
in Appendix B. The latter, denoted <Pa3>, results from taking the orienta- 
tion of the molecules on the underlying FCC sites to be random. The 
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remarkable insensitivity of the Madelung energy (indeed, the total energy) 
to the orientation of the molecules for small values of p/a will be dis- 
cussed below. Further details are given in Appendix A. 

V: Metallic Hydrogen 

The total static Hartree-Fock energy in Ryd/electron for FCC monatomic 

hydrogen using the Abrikosov wavefunction is shown as a function of the 

electron density parameter r^ in Figure 3. Also displayed in this figure 

are the constituent energies--the non-uniform parts of the kinetic, 

electron-proton, electron-electron and exchange energies. For comparison 

a curve showing the homogeneous interacting electron gas Hartree-Fock 

teriTis and the Wigner-Seitz sphere electron-proton energy (i.e., the analog 

of the Madelung energy), is also displayed. Perhaps the most striking 

feature of the H-F curve is how little it departs from this simple estimate: 

the value of r^ at the minimum is essentially unchanged by inclusion of the 

non-uniform terms (the "band structure energy"). It is also noteworthy 

that the two terms in principle hardest to compute--the electron-electron 

and exchange energies--are very small (and their sum even more so) over 

most of the region of interest. An identical calculation using the 

Abrikosov wavefunction for FCC monatomic hydrogen was performed, but over a 

limited range of r^, by Ross and McMahan. We will be interested, however, 

in a discussion of screening over a large density range. At very low density 

28 

(large r^) we find 

ll''i ^el ~ ^ 

s 

where the terms are, respectively, the kinetic, electron-proton. 
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elec^^on-electron , and exchange energies. The mininuni, for AaQ = 27/32 , 
gives a Hartree-Fock energy 

2 

= -0.7119 Ryd/electron (27) 

s 

to wh'^ch the total energy in Figure 3 tends asymptotical ly (and rather 
quickly). As mentioned above, the fact that this value is not that of an 
isolated hydrogen atom reflects the restricted Hartree-Fock substitution 
of a pseudo-atom (an equal mixture of spin-up and spin-down states) for 
the genuine spin-polarized atom, and the total absence of correlation 
energy. As is typ'ical of variational calculations, the energy approaches 
its asymptotic variational minimum faster than do ■'ts variational parameters, 
here only Aa^. 

The inclusion of the correlation energy in the electron gas portion of 

29 

the energy via, for example, the Nozieres-Pines interpolation formula, 
will barely bind FCC metallic hydrogen; the Hartree-Fock minimun energy, 
at r^ = 1.67 and AaQ = 1.24, is -0.929 Ryd/electron: the inclusion of this 

correlation estimate brings the total energy to -1.028 Ryd/electron. The 

H-F energy thus compares very well with other recent H-F calculations for 

30-32 31 

metallic hydrogen, for example, that of Tua and Mahan,. who found 

(using a similar basis of functions but with three variational parameters) 

30 

a H-F energy of -0.9327 Ryd/electron. In addition, Harris et al. found 
a minimizing value of -0.931 Ryd/electron with Aag = 1.2b, using a more 
complicated wavefunction. Ross and McMahan naturally found results 
identical to ours in the density range where comparison can be made. 

Of considerable interest is the density dependence of the wavefunction 
decay parameter A, shown in Figure 4; the original result of Abrikosov, 
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and the Thomas -Fermi prediction 



( 28 ) 


are also displayed for comparison. We may interpret this figure as an indica 
tion of the metallic nature of the monatomic solid as follows. If we con- 
sider the response of the electrons in metallic hydrogen to the fields of 
the protons, then in the Thomas-Fermi model of metallic screening (valid 
as r^-rO) we would expect the induced electron charge density to be 


*^ind'^^ 47rr 


-k. 


TF^ 


(29) 


On the other hand, around an isolated hydrogen atom the charge density is 


p"(r) 


-e 


r 

TT 


■2Xr 


(30) 


with X = l/fiQ- In a condensed metallic phase we would expect X and kjp to 
scale similarly as a function of the density r^. What is remarkable, 
hov ’ver, is the range of density over which the H-F and Thomas-Fermi screen- 
ing parameters actually do track one another; the charge density in this 
model of monatomic hydrogen bears little resemblance to that of an isolated 
atom, and we conclude ^t'at the system is genuinely metallic. In this con- 
text Abrikosov's original calculation is inappropriate for densities of 

33 

physical interest. 


VI: Pa3 Molecular Hydrogen 

In Figure 5 are shown the curves for Pa3 molecular hydrogen analogous 
to those in Figure 3 for the monatomic phase. Here, once again, we display 
a comparison curve, consisting of the interacting electron gas energy and 
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the minimized Madelung energy for the <Pa3> structure, Eq. B5. While there 
is a superficial resemblance to their metallic counterparts, all the 
constituent energies show a remarkable transition at r^ ^2.8, although 
the total energy is quite insensitive to these rather precipitous changes. 

The density dependence of Aag (Figure 6) again closely follows the 
Thomas-Fermi wavevector, as expected for the spherical Fermi surface 
resulting from the use of wavefunction (1). Changes in near r^ = 2.8 

are relatively slight, but a marked flattening-out of beyond r^ 

2.8 is to be noted. 

The density dependence of the interproton spacing (here D/a, i.e., 
the ratio of half the molecular bond length to the lattice constant of the 
conventional simple cubic cell) is shown in Figure 7, and it is clear that 
the abrupt drop in D/a for r^ = 2. 6-3.0 is the origin of the pronounced 
changes in the constituent energies in Figure 5. Shown in the same figure 
are (1) the curves expected if the molecular bond length is frozen at its 
zero pressure value, (ii) the results of the recent density functional 
calculation of Chakravarty et al.,^^ (iii) the results of Liberman, who 
used a modified density functional procedure with a sphericalized potential, 

or 

and (iv) the results of Ramaker, Kumar and Harris for simple cubic oriented 
molecular hydrogen. 

The computed dependence of D/a on density may be understood in terms 
of a qualitative change in the nature of the solid molecular phase around 
r^ = 2.8. For extremely high densities (small r^ ) Figure 2 shows that to 
minimize the Madelung energy, D/a should asymptotical 1y approach a constant 
value ^ 0.27, as we observe in Figure 7. At extremely low densities 
on the other hand (pg-^^), one should recover in an exact treatment a 
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lattice of non-interacting hydrogen molecules with fixed 2D/aQ = 1.401- For 
Pa3 structure this would imply D/a 0.2]7/r^. 

The static Hartree-Fock calculation clearly does rather poorly 
in giving the expected quatitative behavior for low density. However, 
the rather sudden change in behavior does clearly reflect a transition from 
a solid whose essential structural units are "locked into" the unit cell 
as a whole (so that all intracell rlar distances scale together, correspond- 
ing to D/a = constant) to one whose structural units are the molecules 
themselves, with fixed bond length (i.e.. D/a = const/r 2 ). The roughly 
flat dependence of A on r^ beyond r^ % 2.8 in Figure 6 is also charac- 
teristic of weakly interacting molecules. This qualitative change is present, 

in an attenuated form (at a much smaller value of r^) in simple second-order 
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structural expansion results. 

The calculations of Liberman^^ and Chakravarty et a1.^^ use the 
Wigner-Seitz method of sphericalizing the crystalline cell containing one 
molecule and, as such, discard detailed information about the lattice to 
which the total energy is relatively insensitive (though structural parameters 
like D/a may depend sensitively on the approximation). For this reason, 
neither of these calculations exhibit any feature which may be identified as 
the qualitative change we observe. The calculation of Ramaker et al.^® does 
keep a crystalline cell and predicts (at a much higher density) a rapid 
plunge of D/a, but it is for a structure which has not been observed for 
hydrogen. 

It is reasonable to associate the 'transition' at r 2.8 with the 

■ s 
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"rotational transition" marking the onset of free rotation of the hydrogen 
molecules about their centers of mass (ignoring, as usual, the distinction 
between para- and ortho- hydrogen) from the hindered angular oscillation or 
libration characteristic of higher pressures. It is clear that the closer 
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the two protons of a given molecule approach the more nearly spherical is the 
associated charge cloud, so the less important are the configurations of 
other such molecules in ^he same (or other) unit cells. Hence merely on the 
basis of energetics (the insensitivity of the total energy, Figure 5, to 
variation in D/a beyond r^ % 2.S) and proximity (the rapid reduction of the 
bond length shown in Figure 7 around r^ = 2.8), any tendency toward free 
rotation will be markedly increased beyond r^ s 2.8. 

The scaling together of intracellular dimensions below r^ ^ 2.8 is 
typical of metallic systems, where it is the volime-dependent electron gas 
contributions and the interactions of individual ions with the electron gas 
which stabilize the system. In a molecular insulating solid the essential 
entities (the molecules) are neutral and it is fluctuating dipole or higher 
multipolar interactions which stabilize the solid. Simple arguments of Mott 
and others^^’^*^ indicate that for a system with long-range Coulomb interac- 
tions a metal-insulator transition will occur; the system will be metallic 
for an average nunber density n if 

• ( 31 ) 

for which the transition density is r^ n- 3.1. Using the simple cubic 
monatomic hydrogen lattice as a model for such a transition. Rose, Shore 
and Sander^^ do indeed find a metal -i ns ulator transition (from a paramag- 
netic metal for higher density to a ferromagnetic insulator at low density) 
at r^ = 2.84. The precise relevance of results for the monatomic phase to 
the diatomic phase is unclear, but it should be remembered that both systems 
are described by the same Hamiltonian, with an instability toward molecular 
pairing as a structural symmetry breaking. It is believed that as a function 
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of increasing pressure the diatomic phase will first become conducting by 
(indirect) band overlap, and only at higher pressure ^ will the diatwnic 
order disappear. The actual onset of free rotation in (para) hydrogen 

may also be assisted by a Peierls distortion. 

In this context it is important to note that for molecular hydrogen 
the inclusion of lattice dynamics is crucial to obtain a quantitative under- 
standing of the equation of state for low pressures. This problem occurs 
for the monatomic case as well, where the preference of static structural 
expansion calculations for anisotropic planar structures is no longer present 
when a careful self-consistent phonon calculation is made.^ It may be under- 
stood qualitatively by observing that we expect phonon energies per electron 
in the condensed phase to scale roughly as the ionic plasma frequency 
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for m and M the electron and proton masses, respecti v'ily. This quantity 
ranges from about 0.04 Ryd at r = 1 to 0.005 Ryd at r = 4. For compari- 

O W 

son, the zero pressure binding energy per electron for solid molecular 
hydrogen is only about 0.0006 Ryd. The importance of a unified treatment 
of phonons together with the electronic system has been stressed by Chak- 
ravarty et al. For the molecular phase they used a Wigner-Seitz approxi- 
mation; within the W-S sphere around a given molecule they spherical ized 
the ionic potential by smearing the protons out on a shell of radius D, 
reintroducing deviations from spherical symmetry of the total potential within 

the spherical cell by perturbation theory. They found good agreement for molecular 
hydrogen at low density for the binding energy, interproton spacing and 
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optic mode ("stretch'') frequency, The phonon problem manifests itself, 
however, in the fact that the computed equation of state E(rg) showed an 
overall shallow minimifn at r^ =2,1, even when an estimate of the effects 
of phonons was included; the experimental P = 0 density for diatomic 
hydrogen corresponds to r^ = 3,10. 

We may gain a qualitative feeling for the importance of phonons in 
the present calculations of the molecular phase by examining Figure 8, which 
shows (for fixed overall density r^ = 2.8) the dependence of the second 
order structural expansion total energy for the Pa3 and rotational ly averaged 
<Pa3> structures (see Appendix B for a brief description of the second 
order structural expansion and the rotationally averaged <Pa3> structure). 
Although the detailed values of D/a and energy are not reliable, we note 
that (i) the Pa3 curve is extraordi nari ly flat over a very large range of 
D/a out to the simple cubic value D/a = /3/4 (using Hubbard-Geldart-Vosko 
screen! ng--see Appendix B); (ii) the rotationally averaged <Pa3> curve 
falls belo w the Pa3 curve, indicating a possible preference for a lattice 
of free rotators. Moreover, both the splitting between the two curves 
and the "flatness" of the Pa3 curve are of the order of a typical phonon 
energy, indicating that only the inclusion of phonons can specify 
reliably the density-dependence of the molecular inter-proton spacing, the 
onset of the putative rotational solid, and ultimately the low-pressure 
molecular equation of state. 

VII: Sunmary and Conclusions 

The main emphasis . here has been on the Pa3 molecular hydrogen structure. 
For this case screening is found to be quite metallic-like below Ri 2.8; 
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for larger both the screening parameter and the molecular 'bond length' 
cease to depend strongly on density, suggesting a transition to a solid of 
weakly interacting molecules. The detailed nature of this transition has 
yet to be established but a rotational transition and a metal-insulator 
transition at r^ ^ 2,8 are possible candidates. In making the connection, 
it is important to include the effect of phonons. 

Our principal conclusion is that cellular methods which discard 
information about the details of the lattice may do very well when used to 
compute total energies, at the expense, however, of possibly missing 
important structural changes to which the total energy is relatively 
insensitive. We have confirmed once more that it is the difficulty of 
obtaining a reliable first principles molecular equation of state which 
remains the primary obstacle to the understanding of hydrogen at high 
pressure. 
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Ap pendix A : Computational Details 

In the FCC metallic phase (cx^ - “1.7917472304) the only variational 
parameter is the quantity X characterizing the shape of the hydrogen wave- 
function. In practice for each value of r^ eight values of X3 q straddling 
the minimuD (as a function of XaQ) were scanned, a polynomial was fit to 
the resulting energy values, and the appropriate zero of the derivative 
was used as the physically relevant value of Xa^. The energy was then 
evaluated using the curve fit. 

The computational problem for the Pa3 molecular hydrogen phase is of 
a different order of magnitude than for the monatomic phase; not only are 
there two variational parameters, but the lattice itself is of a symmetry 
low enough that systematic exploitation of symmetry is difficult. In order 
to scan a reasonable range of X and D/a in permissible computer time for 
each value of r^ , five values of X3 q and eight of D/a were selected to 
encompass the energy minimum. The resulting table of numbers (from k-space 
sLmmations extending only far out enough to reflect the essential dependence 
on the variational parameters) was interpolated using a bicubic spline 
method, and the minimmi found by a two-dimensional Newton-like method. The 
electronic energy was then re-evaluated at the correspond! ng values of Xag 
and D/a, summing out as far as possible in k-space to improve the accuracy 
of the total energies. 

For several large values of where the total energy is a more 
rapidly varying function of the parameters D/a and XaQ, computer runs with 
a mesh about twice as fine were made to give more' reliable values for the 
structural parameters. ' Even though we feel the variational parameters to 
be accurate (within the limitations of the model), the absolute energies 
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for the molecular phase still contain small uncertainties that make compari- 
son of the equations of state of the molecular and metallic hydrogen diffi- 
cult, especially if the intent is to extract, for example, the transition 
pressure from one phase to the other. Indeed, only the availability of the 
Floating Point Systems Array Processor (and its associated Cornell Fortran 
compiler) made the Pa3 calculation feasible at all in its present form, 
Fluctuations evident in the computed points are due to the sensitive 
dependence of the structural energy on the parameter D/a, and associated 
uncertainty in the output of the fitting and minimization programs. For 
several values of large r^ finer- grid computer runs, indicated by triangles, 
were made. As usual for variational calculations, 'noise' of order 6 in the 
variational parameters emerges as noise in the energy of order 6 . The 
error bars in Figure 6, for example, are indicative of the worst variation 
in results from normal-length computer runs straddling the minimtin in 
parameter space in different ways. 


Appendix B : The Structural Expansion and the Rotational l y Averaged <Pa3> 

Structure 

',‘lithin the structural expansion formalism one begins with a lattice 
of positive ions in the presence of a unifonn compensating negative back- 
ground and, in the same volurie, a uniform interacting electron gas of 
identic?! overall charge density, together with its uniform compensating 
positive background (so that each system is separately neutral). One 
then introduces the electron-ion interaction using perturbation theory. 

To lowest ot'der in the electron-ion interaction one finds 
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where is the energy of the uniforiri interacting electron gas at the 

appropriate density^ is the Madelung energy, given by 
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(where p+^ is the q-th Fourier coinponent of the ion density), and the 
"second order band structure energy" e 1?^ is 
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where v^. is the electron-ion interaction potential (or pseudopotential) 
with e(q) the dielectric function of the uniform interacting electron gas. 
For our calculations we used the Hubbard dielectric function as modified 
by Geldart and Vosko,^^ and the bare Coulomb potential appropriate to 
hydrogen. 

In both expressions above the quantity of interest is the combination 
pi p]->. In the context of an FCC lattice of freely rotating hydrogen 
molecules it is useful to consider an FCC lattice for which the molecular 
bond orientation is a random function of position. One must distinguish 
between the contribution to pi p]^ from the two protons associated with a 
given FCC site (which are perfectly correlated in the molecule there) and 
those associated with other (randomly oriented) molecules. We find 
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where the ^ are the reciprocal lattice vectors for the real space FCC 
lattice; the expression above readily divides into discrete reciprocal 
lattice contributions and a continuous part, which will give rise to integrals 
in the Madelung and band structure energies. Here N is the number of FCC sites, 
Provided the "charged spherical shells" resulting from the rotational 
averaging process do not overlap (for the FCC case this is value for 0 < 

D/a < /?/4) the Madelung energy can then be written 
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where is the Madelung constant of the FCC structure and z = 
The band structure energy, to second order, becomes 
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These expressions were used in the evaluation of the <Pa3> curves (Figures 2 and 8) 
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Figure Captions 

Fig. 1. The a-N 2 or Pa3 structure assuned for the calculation of diatomically 
ordered dense hydrogen. 

Fig. 2. The Madelung constant as a function of interproton spacing (D) (as 
a fraction of conventional cubic cell dimension a). Full line: 

Pa3 structure. Dotted line: rotationally averaged Pa3 (see text). 

Fig. 3. The total stUic Hartree-Fock energy (per electron) for monatomic 
hydrogen obtained with the Abrikosov trial state (1) (lower curve). 
The component parts of the energy are discussed in the text. 

Fig. 4. Density dependence of the variational orbital range parameter X 

in the Abrikosov state (1) for FCC monatomic hydrogen (solid curve). 
The dashed line is the corresponding range in Thomas -Fermi linear 
screening model. The dotted curve is Abrikosov's original cal- 
culation. The numerical details are found in Appendix A. / 

Fig. 5. The total static Hartree-Fock energy (per electron) for diatomically 
ordered hydrogen, minimized with respect to D/a and X. The component 
parts of the energy and their rapid changes near r^ = 2.8 are dis- 
cussed in the text. Note that the total energy is relatively 
featureless at r^ 'v 2.8. 


Fig. 6. Density dependence of the variational orbital range parameter X for 
diatomic hydrogen (solid line). Again the dashed line gives the 
linear-screening result in the Thomas-Fermi approximation (see Fig. 
4). 
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Fig. 7. Minimizing vi>iues of D/a for the diatomically ordered phase of 

hydrogen (solid line and points). Also shown is the curve expected 
if the interproton spacing is held constant at its zero pressure 
value (D = (P = 0) value), the results of Chakravarty et al. 

(Ref. 14), Liberman (Ref. 35), and Ramaker et al. (Ref. 17). 

Fig. 8. Second order structural expansion energies (at r^ = 2.8) for a 

Pa3 structure (solid line) and rotationally averaged Pa3 structure 
(dashed line). Note the extreme insensitivity to D/a (see text). 


